UNIT WTStdHdr;
{***************************************************************}
{*                                                             *}
{*   Copyright (c) Quinn-Curtis, 1992                          *}
{*                                                             *}
{*   Filename:  RTSTDHDR.PAS                                   *}
{*   Revision:  3.0                                            *}
{*   Date:      2/3/92                                        *}
{*                                                             *}
{*   Description: -Routines for virtual arrays                 *}
{*                -Standard numeric types, arrays structures   *}
{*                -Defines BGI Directory used in demos.        *}
{*                                                             *}
{***************************************************************}
INTERFACE

uses  Objects,
      Cartocmds, tit_win;

CONST defaultbgidir = 'c:\TP\BGI';  {defines bgi & chr directory for demos}

TYPE
    REALTYPE = DOUBLE;              {user defined type}
    RealPntr = ^RealType;           {pointer to realtype}
    matstruct = record              {structure for virtual arrays}
                  mr, mc, modef: WORD;
                  mpntr: ^RealPntr;
              end;

{** Available exportable functions **}

FUNCTION GetMatError: INTEGER;
FUNCTION  PowerCalc( realnum: REALTYPE; power :  REALTYPE ):  REALTYPE;
FUNCTION  HeapOverflow(size: longint): BOOLEAN;


TYPE
  P_Byte = ^T_Byte;
  T_Byte = object(TObject)
    data : BYTE;
    constructor Init(NewByte: BYTE);
    destructor Done; virtual;
  end;
  P_Int = ^T_Int;
  T_Int = object(TObject)
    data : Integer;
    constructor Init(NewInteger: Integer);
    destructor Done; virtual;
  end;
  P_Real = ^T_Real;
  T_Real = object(TObject)
    data : REALTYPE;
    constructor Init(NewReal: REALTYPE);
    destructor Done; virtual;
  end;

  POINTER_C = ^TPOINTER_C;
  TPOINTER_C = object(TCollection)

  constructor Init(Size:Integer);
  destructor Done; virtual;

  FUNCTION  GetPntrReal(n: WORD): REALTYPE;
  PROCEDURE PutPntrReal(r:REALTYPE);
  FUNCTION  GetPntrInt(n: WORD): INTEGER;
  PROCEDURE PutPntrInt(i:INTEGER);
  FUNCTION  GetPntrByte(n: WORD): BYTE;
  PROCEDURE PutPntrByte(value: BYTE);
  PROCEDURE Put_PntrReal(n:word; r:REALTYPE);
  PROCEDURE Put_PntrInt(n:word; i:INTEGER);
  PROCEDURE Put_PntrByte(n:word; value: BYTE);

end;

  PROCEDURE DefMat(VAR m: matstruct; r,c: WORD; VAR err: integer);
  PROCEDURE FreeMat(VAR m: matstruct);
  PROCEDURE WrMat(m: matstruct; r,c: WORD; v: REALTYPE);
  FUNCTION  RdMat(m: matstruct; r,c: WORD): REALTYPE;
  PROCEDURE CpyMat2Pntr(m: matstruct; c: word;
                      p: POINTER; n: word);

procedure prod_mat(var m1,m2: matstruct; l,m,n: integer; var m3: matstruct);
procedure scaler_prod_mat(var m1:matstruct; sval:double; numrow,numcol:
                              integer; var m2:matstruct);
procedure add_mat(var m1,m2:matstruct; numrow,numcol:integer; var m3:matstruct);
procedure transpose_mat (var m1:matstruct; numrow,numcol:integer; var m2:matstruct);
function deter_mat(var m1:matstruct; numrow:integer): double;
function invert_mat(var matdata:    matstruct;
                         numcol:     integer;
                     VAR det:        double;
                     VAR invary:     matstruct):integer;

IMPLEMENTATION

VAR sizerealtype, sizereal, sizesingle, sizedouble, sizeinteger,
    sizepointer, sizelong: WORD;
    virtualerr : INTEGER;      {error if non zero }

FUNCTION GetMatError: INTEGER;
BEGIN
  GetMatError := virtualerr;
END;

{raise realnum to power}
FUNCTION PowerCalc( realnum: REALTYPE;
                    power :  REALTYPE ):  REALTYPE;
BEGIN
  PowerCalc := exp(power * ln(abs(realnum)+1.0E-20));
END;

{determines if heap space is available for an allocation of size}
FUNCTION HeapOverflow(size: longint): BOOLEAN;
BEGIN
  IF (size+16) > (MaxAvail) THEN HeapOverflow := TRUE
  ELSE HeapOverflow := FALSE;
END;

{** COLLECTIONS - Ex_POINTER FUNCTIONS **}

(*  T_Byte = object(TObject)   *)
    constructor T_Byte.Init(NewByte: BYTE);
    begin data := NewByte; end;
    destructor T_Byte.Done;
    begin TObject.Done end;

(*  T_Int = object(TObject)  *)
    constructor T_Int.Init(NewInteger: Integer);
    begin data := NewInteger; end;
    destructor T_Int.Done;
    begin TObject.Done; end;

(*  T_Real = object(TObject)  *)
    constructor T_Real.Init(NewReal: REALTYPE);
    begin data := NewReal; end;
    destructor T_Real.Done; 
    begin TObject.Done; end;

(* TPOINTER_C Section *)

constructor TPOINTER_C.Init(Size:integer);
begin TCollection.Init(Size, 5); end;

destructor TPOINTER_C.Done;
begin TCollection.Done; end;

{reads value from element pointed to by p+n}
FUNCTION TPOINTER_C.GetPntrReal(n: WORD): REALTYPE;
VAR index : integer;
BEGIN
  index := n;
  GetPntrReal := REALTYPE(At(index)^);
END;

{write value to element pointed to by p+n}
PROCEDURE TPOINTER_C.PutPntrReal(r:REALTYPE);
BEGIN
  Insert(New(P_Real, Init(r)));
END;

{reads integer value from element pointed to by p+n}
FUNCTION TPOINTER_C.GetPntrInt(n: WORD): INTEGER;
VAR index: INTEGER;
BEGIN
  index := n;
  GetPntrInt := INTEGER(At(index)^);
END;

{write integer value to element pointed to by p+n}
PROCEDURE TPOINTER_C.PutPntrInt(i:INTEGER);
BEGIN
  Insert(New(P_Int, Init(i)));
END;

{reads byte value from element pointed to by p+n}
FUNCTION TPOINTER_C.GetPntrByte(n: WORD): BYTE;
VAR index: integer;
BEGIN
  index := n;
  GetPntrByte := BYTE(At(index)^);
END;

{write byte value to element pointed to by p+n}
PROCEDURE TPOINTER_C.PutPntrByte(value: BYTE);
BEGIN
  Insert(New(P_Byte, Init(value)));
END;

PROCEDURE TPOINTER_C.Put_PntrReal(n:word; r:REALTYPE);
var data: P_Real;
    index:integer;
BEGIN
  index := n;
  data := New(P_Real, Init(r));
  AtPut(index, data);
  Dispose(data, Done);
END;

PROCEDURE TPOINTER_C.Put_PntrInt(n:word; i:INTEGER);
var data: P_Int;
    index:integer;
BEGIN
  index := n;
  data := New(P_Int, Init(i));
  AtPut(index, data);
  Dispose(data, Done);
END;

PROCEDURE TPOINTER_C.Put_PntrByte(n:word; value: BYTE);
var data: P_Byte;
    index:integer;
BEGIN
  index := n;
  data := New(P_Byte, Init(value));
  AtPut(index, data);
  Dispose(data, Done);
END;

{** VIRTUAL ARRAY ROUTINES **}

{** POINTER FUNCTIONS **}
FUNCTION  PntrAddNum (p: POINTER; n, len: INTEGER): POINTER;
TYPE
 LONG = record
   Lp, Hp: WORD;
 END;
VAR memofs, memseg, offset, sum: WORD;
    _result: POINTER;
BEGIN
  offset := (n*len);
  memofs := LONGINT(p) AND $FFFF;
  memseg := LONGINT(p) shr 16;
  sum :=memofs + offset;
  _result := PTR(memseg,sum);
  PntrAddNum := _result;
END;

{reads value from element pointed to by p+n}
FUNCTION GetPntrReal(p: POINTER; n: WORD): REALTYPE;
VAR rpntr: ^REALTYPE;
BEGIN
  GetPntrReal := REALTYPE(PntrAddNum(p,n,sizerealtype)^);
END;

{write value to element pointed to by p+n}
PROCEDURE PutPntrReal(p: POINTER; n: WORD; r:REALTYPE);
VAR rpntr: ^REALTYPE;
BEGIN
  rpntr := PntrAddNum(p,n,sizerealtype);
  rpntr^ := r;
END;

{create virtual array on heap, initial all values to zero}
PROCEDURE DefMat(VAR m: matstruct; r,c: WORD; VAR err: integer);
VAR i,j: WORD;
   rpntr: RealPntr;
   rpntrpntr: ^RealPntr;
   totalsize, columnsize: longint;
BEGIN
  err := 0;
  m.modef := 0;
  IF (r > 4196) AND (c=1) THEN
  BEGIN
    m.modef := 1;
    c := (r div 4196) + 1;
    r := 4196;
  END
  ELSE m.modef := 0;
  columnsize := r*sizerealtype;
  totalsize := columnsize * c + c*4 + 16;
  IF totalsize > MemAvail THEN
    err := 1                             {not enough free heap blocks}
  ELSE
    IF (r>0) AND (c>0) THEN
    BEGIN
      IF HeapOverflow(c*4) THEN
        err := 2                         {not enough contig. heap}
      ELSE
        GetMem(m.mpntr,c*4);     {get heap block of c*4, start addr m.mpntr}
      IF err = 0 THEN
      BEGIN
        m.mr := r;               {matstruct m needs # rows & cols}
        m.mc := c;
        i := 0;
        WHILE ((i < m.mc) and (err = 0)) do
        BEGIN
          IF HeapOverflow(columnsize) THEN
           err := 2
          ELSE
          BEGIN
            GetMem(rpntr,m.mr*sizerealtype);     {get a row of heap}
            rpntrpntr := PntrAddNum(m.mpntr,i,4);
            rpntrpntr^ := rpntr;
            FillChar(rpntr^,m.mr*sizerealtype,0); {initialize to zero}
            INC(i,1);
          END;
        END;
        IF  err <> 0 THEN    {if error occurred, free all heap m was given}
        BEGIN
          FOR j := 0 TO i-1 do
          BEGIN
            rpntrpntr := PntrAddNum(m.mpntr,i,4);
            rpntr := rpntrpntr^;
            FreeMem(rpntr,m.mr*sizerealtype);
          END;
          FreeMem(m.mpntr,m.mc*4);
        END
      END;
    END;
  IF err <> 0 THEN
    _writeln('Insufficient Heap Space for Matrix');
END;

{free memory from heap associated with m}
PROCEDURE FreeMat(VAR m: matstruct);
VAR i: WORD;
   rpntrpntr: ^RealPntr;
   rpntr: RealPntr;
BEGIN
  IF ((m.mc > 0) AND (m.mr > 0))  THEN
  BEGIN
    FOR i := 0 to m.mc-1 do
    BEGIN
      rpntrpntr := PntrAddNum(m.mpntr,i,4);
      rpntr := rpntrpntr^;
      IF rpntr <> NIL
        THEN FreeMem(rpntr,m.mr*sizerealtype);
    END;
    IF m.mpntr <> NIL THEN FreeMem(m.mpntr,m.mc*4);
  END
  ELSE
    virtualerr := 3;
  m.mc := 0;
  m.mr := 0;
END;

{write v to virtual array m at r,c}
PROCEDURE WrMat(m: matstruct; r,c: WORD; v: REALTYPE);
VAR rpntrpntr: ^RealPntr;
    rpntr: RealPntr;
BEGIN
  IF (m.modef = 1) THEN
  BEGIN
    c := r div 4196;
    r := r mod 4196;
  END;
  IF (c-1 < m.mc) AND (r-1 < m.mr) THEN
  BEGIN
    rpntrpntr := PntrAddNum(m.mpntr,c-1,4);
    rpntr := rpntrpntr^;
    PutPntrReal(rpntr,r-1,v);
  END
  ELSE
    virtualerr := 2;
END;

{read value from virtual array m at r,c}
FUNCTION RdMat(m: matstruct; r,c: WORD): REALTYPE;
VAR rpntrpntr: ^RealPntr;
    rpntr: RealPntr;
    _result: Realtype;
BEGIN
  result := 0.0;
  IF (m.modef = 1) THEN
  BEGIN
    c := r div 4196;
    r := r mod 4196;
  END;
  IF (c-1 < m.mc) AND (r-1 < m.mr) THEN
  BEGIN
    rpntrpntr := PntrAddNum(m.mpntr,c-1,4);  {goto column c of m}
    rpntr := rpntrpntr^;
    _result := GetPntrReal(rpntr,r-1);        {get val at rpntr(r,c)}
  END
  ELSE
    virtualerr := 2;
  RdMat := _result;
END;

{write to all values in column c of virtual array m each value in p+i}
PROCEDURE CpyPntr2MatCol(p: POINTER; n: INTEGER;
                         m: matstruct; c: WORD);
VAR i: WORD;
BEGIN
   for i := 0 to n-1 do
     WrMat(m,i,c,GetPntrReal(p,i));
END;

{write to all values in row r of virtual array m
 each value in p+i}
PROCEDURE CpyPntr2MatRow(p: POINTER; n: INTEGER;
                         m: matstruct; r: WORD);
VAR i: WORD;
BEGIN
   for i := 0 to n-1 do
     WrMat(m,r,i,GetPntrReal(p,i));
END;

PROCEDURE CpyMat2Pntr(m: matstruct; c: word;
                      p: POINTER; n: word);
VAR i: WORD;
BEGIN
   for i := 0 to n-1 do
     PutPntrReal(p,i,RdMat(m,i,c));
END;

procedure prod_mat(var m1,m2: matstruct; l,m,n: integer; var m3: matstruct);
var
    i,j,k: integer;
    val : REALTYPE;

begin   {procedure prod_mat}
  for i:=1 to l do
    for j:= 1 to n do
    begin
      val := 0.0;
      for k:= 1 to m do
        val := val + RdMat(m1,i,k) * RdMat(m2,k,j);
      WrMat(m3,i,j,val);
    end;
end;

procedure scaler_prod_mat(var m1:matstruct; sval:double; numrow,numcol:
                              integer; var m2:matstruct);
var i,j: integer;
    val : REALTYPE;

begin  {procedure scaler_prod_mat}
  for i:= 1 to numrow  do
    for j:= 1 to numcol do begin
      val := sval * RdMat(m1,i,j);
      WrMat(m2,i,j,val);
    end;
end;

procedure add_mat(var m1,m2:matstruct; numrow,numcol:integer; var m3:matstruct);
var i,j: integer;
    val : REALTYPE;

begin {procedure add_mat}
  for i:= 1 to numrow  do
    for j:= 1 to numcol do begin
      val := RdMat(m1,i,j) + RdMat(m2,i,j);
      WrMat(m3,i,j,val);
    end;
end;

procedure transpose_mat (var m1:matstruct; numrow,numcol:integer; var m2:matstruct);
var i,j: integer;
    val : REALTYPE;

begin {procedure mat_transpose_mat}
  for i:= 1 to numrow  do
    for j:= 1 to numcol do begin
      val := RdMat(m1,i,j);
      WrMat(m2,j,i,val);
    end;
end;

function deter_mat(var m1:matstruct; numrow:integer): double;
var i,j,k,lx,f,nxt,ier:                      integer;
    pivot,rac,bignum,temp1,temp2,eval :   double;
    val,val1,val2 : REALTYPE;
    DUM : matstruct;

label 999;

begin  {procedure deter_mat}
  deter_mat := -1;
  f := 1;

  DefMat(DUM, numrow, numrow, ier);
  if ier <> 0 then goto 999;

  for i := 1 to numrow do
    for j := 1 to numrow do begin
      val := RdMat(m1,i,j); Wrmat(DUM,i,j,val);
    end;

  for i:=1 to numrow do
  begin
    bignum := 0;
    for k := i to numrow do
    begin
      eval := abs(RdMat(DUM,k,i));
      if eval - bignum > 0 then
      begin
        bignum := eval;
        lx := k;
      end;
    end;
    if i - lx <> 0 then
      f := -f;
    for j := 1 to numrow  do
    begin
      temp1 := RdMat(DUM,i,j);
      val := RdMat(DUM,lx,j);  WrMat(DUM,i,j,val);
      WrMat(DUM,lx,j,temp1);
    end;
    pivot := RdMat(DUM,i,i);
    nxt := i + 1;
    for j := nxt to numrow do
    begin
      val := RdMat(DUM,j,i); rac := val/pivot;
      for k := i to numrow  do begin
        val1 := RdMat(DUM,j,k); val2 := RdMat(DUM,i,k);
        val := val1 - rac * val2;
        WrMat(DUM,j,k,val);
      end;
    end;
  end;
  temp2 := 1;
  for i := 1 to numrow do
     temp2 := temp2 * RdMat(DUM,i,i);

  deter_mat := temp2 * f;

  999: FreeMat(DUM);
end;

function invert_mat(var matdata:    matstruct;
                         numcol:     integer;
                     VAR det:        double;
                     VAR invary:     matstruct):integer;
VAR
                            (* Limited to 256x256 matrices ! *)
  pivlst:                 ARRAY[1..256,1..2] OF integer;
  pivchk:                 ARRAY[1..256] OF boolean;
  i,j,k,l,lerow,lecol,l1: integer;
  piv,t,leval:            double;
  val, val1, val2 : REALTYPE;

label 999;

PROCEDURE swap(VAR s1, s2: double);
VAR
  temp: double;

BEGIN         {procedure swap}
  temp := s1;
  s1 := s2;
  s2 := temp;
END;

BEGIN {procedure invert_mat}

invert_mat := -1;

det := 1;
FOR i := 1 TO numcol DO
  BEGIN
    pivchk[i] := false;
    FOR j := 1 TO numcol DO begin
      val := RdMat(matdata,i,j);
      WrMat(invary,i,j,val);
    end;
  END;

  FOR i:= 1 TO numcol DO
  BEGIN
    leval := 0.0;
    FOR j := 1 TO numcol DO
       IF NOT(pivchk[j]) THEN
         FOR k := 1 TO numcol DO
           IF NOT(pivchk[k])  THEN
             IF abs(RdMat(invary,j,k)) > leval THEN
             BEGIN
               lerow := j;
               lecol := k;
               leval := abs(RdMat(invary,j,k));
             END;

     pivchk[lecol] := true;
     pivlst[i,1] := lerow;
     pivlst[i,2] := lecol;

     IF lerow <> lecol THEN
     BEGIN
        det := -det;
        FOR l := 1 TO numcol DO begin
          val1 := RdMat(invary,lerow,l);
          val2 := RdMat(invary,lecol,l);
(*          swap(invary^[lerow,l],invary^[lecol,l]);  *)
          WrMat(invary,lerow,l,val2);
          WrMat(invary,lecol,l,val1);
        end;
     END;

     piv := RdMat(invary,lecol,lecol);
     det := det * piv;
     IF det > 1.0e+38 THEN goto 999; (* Prevents overflow *)
     val := 1.0; WrMat(invary,lecol,lecol,val);
     FOR l := 1 TO numcol DO begin
       val := RdMat(invary,lecol,l)/piv;
       WrMat(invary,lecol,l,val);
     end;

     FOR l1 := 1 TO numcol DO
     BEGIN
       IF l1 <> lecol THEN
       BEGIN
         t := RdMat(invary,l1,lecol);
         val := 0.0; WrMat(invary,l1,lecol,val);
         FOR l := 1 TO numcol DO begin
           val := RdMat(invary,l1,l) - RdMat(invary,lecol,l) * t;
           WrMat(invary,l1,l,val);
         end;
       END
     END
   END;

 FOR i := 1 TO numcol DO
 BEGIN
   l := numcol - i + 1;
   IF pivlst[l,1] <> pivlst[l,2] THEN
   BEGIN
     lerow := pivlst[l,1];
     lecol := pivlst[l,2];
     FOR k := 1 TO numcol DO begin
          val1 := RdMat(invary,k,lerow);
          val2 := RdMat(invary,k,lecol);
(*       swap(invary^[k,lerow],invary^[k,lecol])  *)
          WrMat(invary,k, lerow,val2);
          WrMat(invary,l, lecol,val1);
     end;
   END
 END;

 invert_mat := 0;

999:END;

BEGIN
  sizerealtype := sizeof(realtype);
  sizereal     := sizeof(real);
  sizesingle   := sizeof(single);
  sizedouble   := sizeof(double);
  sizeinteger  := sizeof(integer);
  sizepointer  := sizeof(pointer);
  sizelong     := sizeof(longint);
  virtualerr := 0;
END.
